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The global asymptotic nonlinear behavior of some standard iterative procedures in solving 
nonlinear systems of algebraic equations arising from four implicit linear multistep methods 
(LMMs) in discretizing three models of 2 x 2 systems of first-order autonomous nonlinear 
ordinary differential equations (ODEs) is analyzed using the theory of dynamical systems. The 
iterative procedures include simple iteration and full and modified Newton iterations. The 
results are compared with standard Runge-Kutta explicit methods, a noniterative implicit 
procedure, and the Newton method of solving the steady part of the ODEs. Studies showed 
that aside from exhibiting spurious asymptotes, all of the four implicit LMMs can change the 
type and stability of the steady states of the differential equations (DEs). They also exhibit a 
drastic distortion but less shrinkage of the basin of attraction of the true solution than standard 
nonLMM explicit methods. The simple iteration procedure exhibits behavior which is similar 
to standard nonLMM explicit methods except that spurious steady-state numerical solutions 
cannot occur. The numerical basins of attraction of the noniterative implicit procedure mimic 
more closely the basins of attraction of the DEs and are more efficient than the three iterative 
implicit procedures for the four implicit LMMs. Contrary to popular belief, the initial data 
using the Newton method of solving the steady part of the DEs may not have to be close to the 
exact steady state for convergence. These results can be used as an explanation for possible 
causes and cures of slow convergence and nonconvergence of steady-state numerical solutions 
when using an implicit LMM time-dependent approach in computational fluid dynamics. 


1. Background and Objective 

It has been shown recently by the authors and 
others [Yee et al. , 1991; Yee & Sweby, 1992; 
Lafon & Yee, 1991; Lafon & Yee, 1992; Griffiths 
et al. , 1992a, b; Sweby Sz Yee, 1991; Yee et al , 
1992; Mitchell & Griffiths, 1985; Iserles, 1988; 
Iserles et al. , 1990; Stuart, 1989; Dieci & Estep, 


1990] that the dynamics of the numerical dis- 
cretizations of nonlinear differential equations 
(DEs) can differ significantly from that of the 
original DEs themselves. For example, it was shown 
in Yee et al. [1991], Yee & Sweby [1992], and Lafon 
& Yee [1991, 1992] that the discretizations can pos- 
sess spurious steady-state numerical solutions and 
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spurious asymptotes which do not satisfy the orig- 
inal DEs. These spurious numerical solutions may 
be stable or unstable and may occur both below and 
above the linearized stability limit of the numer- 
ical scheme (on the time step for the equilibrium 
or asymptote of the DE). In Yee k Sweby [1992], 
Lafon k Yee [1991], Sweby k Yee [1991], and Yee 
et al. [1992a] we showed how “numerical” basins 
of attraction can complement the bifurcation dia- 
grams in obtaining the global asymptotic behav- 
ior of numerical solutions for nonlinear DEs. We 
showed how in the presence of spurious asymptotes 
the basins of the true stable steady states could be 
altered by the basins of the spurious stable and un- 
stable asymptotes. One major consequence of this 
phenomenon which is not commonly known is that 
this spurious behavior can result in a dramatic dis- 
tortion and, in most cases, a dramatic shrinkage 
and segmentation of the basin of attraction of the 
true solution for finite time steps. Such distortion, 
shrinkage, and segmentation of the numerical basins 
of attraction will occur regardless of the stability of 
the spurious asymptotes. 

We use the term “spurious asymptotic numer- 
ical solutions” to mean asymptotic solutions that 
satisfy the discretized counterparts but do not sat- 
isfy the underlying ordinary differential equations 
(ODEs) or partial differential equations (PDEs). In 
other words, asymptotic solutions that are asymp- 
totes of the discretized counterparts but are not 
asymptotes of the DEs. Asymptotic solutions here 
include steady-state solutions, periodic solutions, 
limit cycles, chaos, and strange attractors. Here, 
the basin of attraction is a domain of a set of ini- 
tial conditions whose solution curves (trajectories) 
all approach the same asymptotic state. Also, we 
use the term “exact” and “numerical” basins of at- 
traction to distinguish “basins of attraction of the 
underlying DEs” and “basins of attraction of the 
discretized counterparts”. 

In view of the above nonlinear behavior of nu- 
merical schemes, it is possible that numerical com- 
putations may converge to an incorrect steady state 
or other asymptotes which appear to be physically 
reasonable. One major implication is that what is 
expected to be physical initial data associated with 
a true steady state might lead to a wrong steady 
state, a spurious asymptote, or a divergence or non- 
convergence of the numerical solution. In addition, 
the existence of spurious limit cycles [Yee et a/., 
1992; Yee et a/., 1991; Yee k Sweby, 1992] may re- 
sult in the type of nonconvergence of steady-state 


numerical solutions observed in time-dependent 
approaches to the steady states. It is our belief 
that the understanding of the symbiotic relation- 
ship between the strong dependence on initial data 
and permissibility of spurious stable and unstable 
asymptotic numerical solutions at the fundamental 
level can guide the tuning of the numerical param- 
eters and the proper and/or efficient usage of nu- 
merical algorithms in a more systematic fashion. It 
can also explain why certain schemes behave nonlin- 
early in one way but not another. Here, strong de- 
pendence on initial data means that for a finite time 
step At that is not sufficiently small, the asymptotic 
numerical solutions and the associated numerical 
basins of attraction depend continuously on the ini- 
tial data. Unlike nonlinear problems, the associated 
numerical basins of attraction of linear problems are 
independent of At as long as At is below a certain 
upper bound. 

Studies in Yee et al. [1991], Yee k Sweby [1992], 
Lafon k Yee [1991, 1992], Griffiths et al. [1992a, b], 
Sweby k Yee [1991], Yee et al [1992], Mitchell k 
Griffiths [1985], Iserles [1988], Iserles et al [1990], 
Stuart [1989], and Dieci k Estep [1990] are par- 
ticularly important for computational fluid dynam- 
ics (CFD), since it is a common practice in CFD 
computations to use a time-dependent approach to 
obtain steady-state numerical solutions of compli- 
cated steady fluid flows which often consist of stiff 
nonlinear PDEs of the mixed type. When a time- 
dependent approach is used to obtain steady-state 
numerical solutions of a fluid flow or a steady PDE, 
a boundary value problem is transformed into an 
initial-boundary value problem with unknown ini- 
tial data. If the steady PDE is strongly nonlinear 
and/or contains stiff nonlinear source terms, phe- 
nomena such as slow convergence, nonconvergence 
or spurious steady-state numerical solutions com- 
monly occur even though the time step is well below 
the linearized stability limit and the initial data are 
physically relevant. 

It is also a common practice in CFD to use im- 
plicit methods to solve stiff problems. Such schemes, 
however, introduce the added difficulty of solving 
the implicit equations (nonlinear algebraic equa- 
tions) in order to obtain the solution at the next 
time level. Various options such as linearization or 
iteration are available for this purpose. In Yee k 
Sweby [1992], we included a study on the dynamics 
of a noniterative linearized version of the implicit 
Euler and trapezoidal methods. In this work we 
generalize our earlier work [Yee k Sweby, 1992] to 
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include iterative solution procedures of the implicit 
discretized equations, namely, simple iteration and 
full and modified Newton iteration. In addition, we 
also study the 3-level backward differentiation for- 
mula (BDF) and a midpoint implicit method (one- 
legged method of the trapezoidal method) with 
these four methods (noniterative versus iterative) of 
solving the resulting nonlinear algebraic equations 
applied to it. A comparison of the above combina- 
tion of implicit methods with the nine explicit meth- 
ods studied in Yee Sz Sweby [1992], and with the 
“straight” Newton method in solving the steady 
part of the equations, is also performed. We use 
the term “straight Newton” to distinguish it from 
the combination of “implicit LMM -I- Newton” type 
of solution procedure. Some study on variable time 
step control to avoid spurious asymptotes will also 
be addressed. A condensed version of this work 
will appear in Yee &: Sweby [1993]. Generaliza- 
tion of the study to include grid adaptation as one 
of the sources of nonlinearity and/or stiffness that 
is introduced into the discretized system will be 
reported in Sweby & Yee [1994] and Budd et al. 
[1994]. 

2. Relevance and Outline 

Since all four implicit methods under consideration 
are linear multistep methods (LMMs) they will not 
exhibit spurious steady states (fixed points of order 
one). However, as discussed in Yee &; Sweby [1992] 
and as we shall see in later sections, some implicit 
LMMs have the property of increasing the stability 
range for the stable fixed points of the ODE [Dieci 
& Estep, 1990], accompanied in some instances by 
the stabilization of unstable fixed points of the ODE 
[Yee & Sweby, 1992]. In addition, the method of so- 
lution of the implicit equations generated by these 
schemes can itself contribute to the dynamics of 
the discretization since different numerical methods 
and/or solution procedures result in entirely differ- 
ent nonlinear discrete maps. Iserles [1988] and Dieci 
& Estep [1990] were the first to examine some of the 
stability issues. Our attempt here is to address is- 
sues that were not investigated in Iserles [1988] and 
Dieci & Estep [1990], and in our companion paper 
[Yee & Sweby, 1992]. Our main purpose is to study 
the global asymptotic behavior in terms of bifurca- 
tion diagrams and numerical basins of attraction of 
these four procedures for solving nonlinear systems 
of algebraic equations arising from implicit LMM 
discretizations. 


Although we purposely selected the model equa- 
tions with known analytical solutions, depending on 
the scheme, the dynamics of their discretized coun- 
terpart are very difficult and might not be possi- 
ble to analyze analytically. Only some analysis is 
possible. Part of the global asymptotic numerical 
solution behavior can be obtained by the pseudo 
arclength continuation method devised by Keller 
[1977], a standard numerical method for obtaining 
bifurcation curves in bifurcation analysis. Besides 
not being able to provide the numerical basins of 
attraction, one deficiency of the pseudo arclength 
continuation method is that, for problems with 
complicated bifurcation patterns, it cannot pro- 
vide the complete bifurcation diagram without a 
known solution for each of the main bifurcation 
branches. For spurious asymptotes it is usually not 
easy to locate even just one solution on each of these 
branches. For the majority of the cases where rig- 
orous analysis is impractical, we utilized numeri- 
cal experiments. Also, analytical representations 
(except in isolated cases) for numerical basins of 
attraction rarely exist for nonlinear DEs. Methods 
such as generalized cell mapping [Hsu et a/., 1982; 
Flashner &; Guttalu, 1988] can provide an efficient 
approach to locating these basins, but might not be 
exact. Here, our aim is to numerically compute the 
basins of attraction as accurately as possible and in 
the most straightforward way in order to illustrate 
the key points. 

Due to the complicated nature of these dis- 
crete maps, analysis without a supercomputer is 
nearly impossible. The nature of our studies re- 
quires the performance of a very large number of 
simulations with different initial data; this can be 
achieved by use of the highly parallel Connection 
Machines (CM-2 or CM-5) whereby each processor 
could represent a single initial datum and thereby 
all the computations can be done in parallel to pro- 
duce detailed global stability behavior and the re- 
sulting basins of attraction. With the aid of highly 
parallel Connection Machines, we were able to de- 
tect a wealth of the detailed nonlinear behavior of 
these schemes which would have been overlooked 
had isolated initial data been chosen on the Cray- 
YMP or other serial or vector machine. 

Outline: 

The outline of this paper is as follows. Section 3 
reviews background material. Section 4 describes 
the three 2x2 systems of nonlinear first-order au- 
tonomous model ODEs. Section 5 describes the 
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four implicit LMM time discretizations and the four 
commonly used noniterative and iterative proce- 
dures in solving the resulting nonlinear algebraic 
equations. Section 6 discusses the combined bifur- 
cation diagrams and numerical basins of attraction 
of the underlying schemes for the three model equa- 
tions. This section also includes the result of using 
straight Newton method to solve the steady part 
of the ODE ^ = S(U ) (i.e., solving for the solu- 
tions of S(U) = 0). Due to the complexity of the 
subject matter most of our study used a fixed time 
step. A brief discussion on variable step size con- 
trol will be included. The paper ends with a sum- 
mary and a discussion of the implications of the 
numerical studies in Sec. 7. Our studies reveal that 
unconditionally stable LMMs perform better than 
standard explicit methods and that noniterative im- 
plicit schemes are more efficient and mimic more 
closely the true behavior of the governing equations 
than the iterative implicit procedure for all of the 
four studied implicit LMMs. Studies further show 
that standard variable time step control, although 
performing better with very restricted time steps 
than the constant time step case, did not allevi- 
ate completely the occurrence of spurious numerical 
asymptotes. 

3. Preliminaries 

Consider a 2 x 2 system cf first-order autonomous 
nonlinear ODEs of the form 

^ =S(U), (3.1) 

where U and S are vector functions of dimension 2, 
and S(U) is nonlinear in U . A fixed point (steady 
state) Ue of an autonomous system (3.1) is a con- 
stant solution of (3.1); that is 

S(U E ) = 0, (3.2) 

where the subscript “E” stands for “exact” and Ue 
denotes the fixed points of the ODE as opposed to 
the additional fixed points of the discretized coun- 
terparts (spurious fixed points) due to the numerical 
methods which we will encounter later. 

Consider a nonlinear discrete map from a finite 
discretization of (3.1) 

U n+1 = U n + D(U n , At), (3.3) 

where At is the time step and D(U n , At) is 
linear or nonlinear in At depending on the nu- 


merical method. A fixed point U of (3.3) is defined 
by U n+l = U n , or 

U = U + D(U, At) (3.4a) 

or 

D(U, At) = 0 . (3.4b) 

A fixed point U of period p > 0 of (3.3) is defined 
by U n+P = U n where U n+k # U n for Jfc < p. In the 
context of discrete systems, the term “fixed point” 
without indicating the period means “fixed point of 
period 1” or the steady-state solution of (3.3). Here, 
we use the term asymptote to mean a fixed point of 
any period, a limit cycle (invariant set), chaos, or a 
strange attractor. 

The type of finite discretization of (3.1) repre- 
sented in (3.3) assumed the use of two-time level 
schemes. Otherwise the vector dimension of (3.3) 
would be 2(k — 1) instead of 2 where k is the num- 
ber of the time level of the scheme. Here, the vec- 
tor function D is assumed to be consistent with the 
ODE (3.1) in the sense that fixed points of the ODE 
are fixed points of the scheme; however, the re- 
verse need not hold. Also, spurious asymptotes that 
are asymptotic numerical solutions of (3.3) but not 
(3.1) can exist, depending on the numerical method 
and A t. It is these features, accompanied by other 
added dynamics, that cause the discretized coun- 
terparts of the underlying ODE to possess a much 
richer dynamical behavior than the original ODE 
which forms the core of this study. Thus, the fixed 
points U of D(U , At) = 0 may be true fixed points 
Ue of (3.1) or spurious fixed points t/5. The spuri- 
ous fixed points Us are not roots of S(U) = 0. That 
is, S(Us) #0. 

Letting U n = U + 6 n , then a perturbation anal- 
ysis on (3.3) by discarding terms of 0(6 2 ) yields 

^._(, + «2§M)*V (3.5) 

Assuming dD (^ At ) ^ q, then the fixed point U is 

stable if the eigenvalues of lie inside the 

unit circle. If both eigenvalues are real and both 
lie inside (outside) the unit circle, then the fixed 
point is a stable (unstable) node. If one is inside 
the unit circle and the other outside, then the fixed 
point is a saddle. If both eigenvalues are complex, 
then the fixed point is a spiral. If the eigenvalues 
lie on the unit circle, then the fixed point of (3.3) 
is indeterminant and additional analysis is required 
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to determine the true behavior of (3.3) around this 
type of fixed point. For a more refined definition 
and the difference in fixed point definition between 
ODEs and discrete maps, see Guckenheimer &; 
Holmes [1983] and Hale &; Kocak [1991] and ref- 
erences cited therein. 

An important feature which can arise (for both 
systems of ODEs (3.1) and their discretizations) 
as the result of a Hopf bifurcation is a limit 
cycle where the trajectory traverses a closed curve 
in phase space. In all but a few simple cases, such 
limit cycles are beyond analysis. It is possible that 
the ODEs posses no limit cycle, but depending 
on the numerical methods, spurious limit cycles 
can be present for the discretized counterparts. 
It is this phenomena that can contribute to the 
nonconvergence of numerical schemes in time- 
dependent approaches to the steady states. 


4. Model Nonlinear First-Order 
Autonomous ODEs 

Three of the four 2x2 systems of nonlinear first- 
order autonomous model ODEs considered in Yee 
& Sweby [1992] are considered here. As before, we 
do not treat any system parameter present in the 
DEs as a bifurcation parameter, but instead keep it 
constant throughout each numerical calculation so 
that only the discretization parameters come into 
play. The systems considered with U T = (u, v ) or 
^ = u -f iv are a 

1. Dissipative complex model: 

^ = z(i + £ - \z\ 2 ) ] (4.1) 


The perturbed Hamiltonian model can be re- 
lated to the numerical solution of the viscous Burg- 
ers’ equation with no source term: 


du 1 d(u 2 ) _ d 2 u 
dt + 2 dx dx 2 


0 > 0 . 


(4.4) 


Let Uj(t) represent an approximation to u(xj , t) 
of (4.4) where xj = j Ax, j = 1 with Ax 
the uniform grid spacing. Consider the three-point 
central difference spatial discretization with peri- 
odic condition uj+j = Uj, and assume Ylj-i u j = 

constant, which implies that %2j= l ~3t~ ~ 0* ^ we 
take J = 3 and Ax = 1/3, then, with e = 9/3 
this system can be reduced to the 2x2 system of 
first-order nonlinear autonomous ODEs (4.3) with 
JJ T = (ui, U 2 ) = (u, v)- In this case, the nonlinear 
convection term is contributing to the nonlinearity 
of the ODE system (4.3). 

These three equations were selected to bring 
out the dynamics of numerics for three different 
types of solution behavior of the nonlinear ODEs. 
The dissipative complex system (4.1) possesses ei- 
ther a unique stable fixed point or limit cycle with 
an unstable fixed point depending on the value of 
e . This is the rare situation where the analyti- 
cal expression of a limit cycle can be found. The 
predator-prey model (4.2) exhibits multiple stable 
fixed points. The perturbed Hamiltonian model 
(4.3), which arises as a gross simplification of the fi- 
nite discretization of the viscous Burgers 1 equation, 
exhibits a unique stable fixed point. Following the 
classification of fixed points of (3.1) in Sec. 3, one 
can easily obtain the following: 


2. Predator-Prey model: 


du 

dt 


= —3 u + 4 u l — 0.5uu — *u J 


dv 

— — -2.lt; T uv ; 
dt 


(4.2) 


3. Perturbed Hamiltonian system model: 

[l - 2 u + u 2 - 2v(l - u) 


du /„ x 3 


(4.3) 


dv n ? \ 3 


j^l - 2v + v 2 - 2u(l - u) . 
Here, e is the system parameter for (4.1) and (4.3). 


Fixed Points of (4-1)' 

The dissipative complex model has a unique fixed 
point at ( u , v) = (0, 0) for e < 0. The fixed point 
is a stable spiral if e < 0 and a center if e = 0. 
For e > 0, the fixed point (0, 0) becomes unstable 
with the birth of a stable limit cycle with radius 
equal to s/e centered at (0,0). Figure 4.1 shows 
the phase portrait ( u-v plane) of system (4.1) for 
e = — 1 and e = 1, respectively. Here, the entire 
(u, v) plane belongs to the basins of attraction of 
the stable fixed point (0,0) if e < 0. On the other 
hand, if e > 0, the entire (u, v) plane except the 
unstable fixed point (0,0) belongs to the basins of 
attraction of the stable limit cycle. 
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Phase Portraits & Basins of Attraction 
Dissipative Complex Equation 




u u 

Fig. 4.1. The phase portrait (u-v plane) of system (4.1) for e = — 1 and e = 1, respectively. 


Phase Portrait & Basins of Attraction 
Predator-Prey Equation 
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U 

Fig. 4.2. The phase portrait and their corresponding basins 
of attraction for system (4.2). 


Fixed Points of (4-%) : 

The predator-prey equation has four fixed points 
(0,0), (1,0), (3,0), and (2.1,1.98). By looking at 
the eigenvalues of the Jacobian of S', one finds that 
(0,0) is a stable node, (2.1,1.98) is a stable spiral, 
and (1,0) and (3,0) are saddles. Figure 4.2 shows 
the phase portrait and their corresponding basins 
of attraction for system (4.2). The different shades 
of grey regions represent the various basins of 
attraction of the respective stable fixed points. 
The white region represents the basin of divergent 
solutions. Note that the trajectories near the un- 
stable separatrices actually do not merge with the 
unstable branch of separatrices, but only appear 
to merge due to the thick drawings of the solution 
trajectories. 

Fixed Points of (4-3): 

The perturbed Hamiltonian (semidiscrete system of 
the viscous Burgers’ equation with three-point cen- 
tral difference in space) has four steady-state solu- 
tions of which three are saddles and one is a sta- 
ble spiral at (1/3, 1/3) for e > 0. For e = 0, the 
stable spiral becomes a center. Figure 4.3 shows 
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Phase Portraits & Basins of Attraction 


Viscous Burger's Equation (Central Difference in Space) 

£ = 0 £ = 0.1 



Fig. 4.3. The phase portrait and their corresponding basins of attraction for system (4.3). 


the phase portrait and their corresponding basins 
of attraction for system (4.3). The shaded region 
represents the basins of attraction for the fixed 
point (1/3, 1/3) for e = 0 and e = 0.1. The white 
region represents the basin of divergent solutions. 
From here on we refer to (4.3) also as a viscous 
Burgers’ equation with central difference in space. 

5. Numerical Methods 

The four LMMs for (3.1) considered here are 

1. Implicit Euler method 

{/”+! = U n + AfS n+1 ; (5.1) 

2. Trapezoidal method 

U n+1 = U n + ^At{S n + S n+1 ) ; (5.2) 

3. 3- Level backward differentiation 
formula (BDF) 

U n+ 1 = U n + -A tS n+1 + ~{U n - U n ~ l ) ; (5.3) 
3 3 


4. Mid-point implicit method (one-legged 
trapezoidal) 

U n+l = U n + Ats ~(U n+1 + U n ) . (5.4) 

Performing standard perturbation analysis on the 
above equations at the fixed points U of the ODE 
system by writing U n = U + S n where S(U) = 0 
and discarding terms of 0(8 2 ) yields 

(5 n+1 = K(U)6 n (5.5) 

where the matrix K(U) is defined implicitly for the 
3-level BDF method. The stability of the pertur- 
bation is therefore governed by the eigenvalues /i of 
the matrix K(U). 

For the implicit Euler method we find that 

K(U) = [I - AtJiU)}- 1 , (5.6) 

where J(U) is the Jacobian dSjdU evaluated at the 
fixed point, while for the trapezoidal and midpoint 
implicit methods 

K(U) =\l- ±AtJ(U)] 1 [/ + l - At J(U)} , (5.7) 
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and for the 3-level BDF method we have the re- 
lationship 

/ - | At J(Z7)] K(U) 2 - 1k(U) + \l = 0 . (5.8) 


Thus, if the eigenvalues of J(U ) are A, we have the 
eigenvalues fi of K(U) as follows: 


Implicit Euler 

i 

1-AtA ’ 

(5.9a) 

Trapezoidal 
Midpoint Implicit 

2 + AtX 
M ~2 - MX ’ 

(5.9b) 

3-level BDF 

1 

** 2 Vl+ 2AtA ’ 

(5.9c) 


with both matrices sharing the same eigenvectors. 
(Note that since the expression (5.8) is a quadratic 
in K(U) there are four possible modes correspond- 
ing to the different eigenvalues A and the two pos- 
sible signs of the square root. However, it is the 
modes taking the negative square root which have 
larger modulus and therefore govern stability of the 
perturbation.) 

The stability of the corresponding fixed points 
based on the eigenvalues of K(U) can be determined 
exactly and were used to check our numerical com- 
putations later. They are summarized in Table 1. 
When numerically computing the full discretized 
equations there are various options which can be 
used to solve Eqs. (5.1), (5.2), or (5.3). We consider 
here linearization (a noniterative procedure [Yee & 
Sweby, 1992]), simple iteration, Newton iteration, 
and modified Newton iteration. 

Linearization is achieved by expanding S n+1 
as S n + J(U n ){U n ^ - t/ n ). Thus, the linearized 
implicit Euler method is 

jyn+l = V n + _ At J{U n )\ ~ l S n , (5.10) 


Straightforward perturbation analysis shows 
that while the matrices K{U) for the perturbed 
form of the implicit scheme differ from the fully 
implicit schemes, their eigenvalues are the same. 
Note, however, that the dynamics away from the 
fixed points need not be identical since the per- 
turbation analysis represents only the local behav- 
ior of hypberbolic fixed points of the unperturbed 
systems. 

Simple Iteration is the process in which, 
given a scheme of the form 

U n+1 =G(U n ,U n+1 ), (5.13) 

we perform the iteration 

Ufc\) = G(U n , U^ 1 ) i/ = l,... (5.14) 

where Uq +1 = U n and U (u) n indicates the iteration 
index. The iteration is continued either until some 
tolerance between iterates is achieved, i.e., 

Il^-^rntol (5.15) 

or a limiting number of iterations, in our case 15, 
have been performed. The major drawback with 
simple iteration is that for guaranteed convergence 
the iteration must be a contraction, i.e., 

II G(U n , V ) - G(U n , W)\\ < a\\V - W\\ , (5.16) 

where a < 1. Whether or not the iteration is a con- 
traction at the fixed points will influence the stabil- 
ity of that fixed point, overriding the stability of the 
implicit scheme. Away from the fixed points the in- 
fluence will be on the basins of attraction. For the 
four LMMs at the fixed points this translates as 
follows: 

Implicit Euler At|| J{U)\\ < a < 1 , (5.17a) 


the Linearized Trapezoidal method (and also the Trapezoidal & -At\\ J(U)\\ <a<l (5.17b) 

Unearned midpoint implicit method) is Midpoint Implicit 2 ~ 


1 -l-i 

U n+1 = U n + At I - -AtJ(U n ) S n , (5.11) 

Li 

and the Linearized 3-level BDF method is 


3-level BDF ^At\\J(U)\\<a<l . (5.17c) 

O 

As we shall see in Secs. 6 and 7, our numerical 
results illustrate this limitation well. 


jjn+l __ jjn _|_ j 


;AtJ(U n ) 


Newton Iteration for the implicit schemes is 
of the form 


x -A tS n + -{U n -U r 

u o 


(5.12) 


rm+l 

U (^+D 


F'(U n , U J+ 1 )- 1 

xF(U n ,U$ l ) 


Note that all three have possible singularities. 
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Table 1. Stability of the fixed points of the perturbed discretized equations. 


System/Fixed Point 

Exact 

Implicit Euler 

Trap 

BDF 

Dissipative Complex 






(0, 0) e < 0 

SS 

SS 


ss 

Stable 

(0, 0) e > 0 

US 

US 

(%?,) 

US 

Unstable (0, ?) 



SS 

(■+.»■“) 


Stable (?, oo) 

Predator Prey 






(0, 0) 

SN 

SN 


SN 

Stable 

(1.0) 

S 

s 

(0, 1) 

S 

Unstable (0, 2) 



SN 

(1, oo) 

S 

Stable (2, oo) 

(3, 0) 

S 

s 

(°4) 

S 

Unstable ^0, 



SN 

(£'“) 


stable (o*°°) 

(2.1, 1.98) 

SS 

ss 


SS 

Stable 

Perturbed Hamiltonian 






G'sH 

SS 

ss 


SS 

Stable 


Key: SN — Stab le Node, SS — Stable Spiral, US — Unstable Spiral, S — Saddle, Trap — trapezoidal, 
7 = f + yje 2 + 4. Stable/Unstable indicates that type could not be determined. Intervals indicate ranges of 
At if not ( 0 , oc). 


where Uq+ 1 — U n . The differentiation is with re- 
spect to the second argument and the scheme has 
been written in the form (for two-time level 
schemes) 

F(U n , t/ n+1 ) = 0 . (5.19) 

Modified Newton Iteration is the same as (5.18) 
except it uses a frozen Jacobian F f {U n , U n ). The 
same tolerance and maximum number of itera- 
tions used for the simple iteration are also used 
for the Newton and modified Newton iterations. In 
all of the computations, the starting scheme for the 
3-level BDF is the linearized implicit Euler. 

We also considered two variable time step con- 
trol methods. The first one is “implicit Euler + 
Newton iteration with Local Truncation Error Con- 
trol” [Dieci & Estep, 1990] 


E?„V„ = VJ*' + [/ - 

X Pm' -V"- At"S(C/,”t 1 )| . = 1,. ■ ■ 

(5.20a) 


with 

Af n = 0.9Ai n_1 


toll 

|| [ /n_[/n-l_A<"-lS({/n) 


(5.20b) 


where the (n + l)th step is rejected if || U n - 
U n ~ l — At” _1 S(f/ n )|| > 2toli- In this case, we set 
A t n ~ l — A t n . The value “toll” is a prescribed tol- 
erance and the norm is an infinity norm. The second 
one is the popular “ode23” method 

ki = S(U n ), 

k 2 = S(U n + At n h), 

k 3 = S{U n + At n (h + fc 2 )/4) , (5.21a) 

U n +1 = u n + At n (k! + k 2 + 4/c 3 )/6 , 

A U n+l = At n (ki +k 2 - 2fc 3 )/3 , 
with 

At "= 09A1 "'‘/ii^ir' (5 - 21b) 
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where the (n + l)th step is rejected if ||A[/ n+1 || > 
toh max{l, ||[/ n+1 ||}. In that case, we set At 71 ' 1 — 
At 71 . Again, “toll” is a prescribed tolerance and the 
norms are infinity norms. 

We also employ “straight” Newton’s method 
in solving for the solutions of S(U) — 0 which is 
the one-step Newton iteration of the implicit Euler 
method of (5.18). 

6. Numerical Results 

As mentioned before, the nature of our calculations 
requires thousands of iterations of the same equa- 
tion with different ranges of initial data on a pres- 
elected (u, v) domain and ranges of the discretized 
parameter space At. The NASA Ames CM-2 and 
CM-5 allow vast numbers (typically 65,536) of cal- 
culations to be performed in parallel and is ideal 
for our studies. Two different representations of the 
numerical basins of attraction are computed on the 
Connection Machines. One is bifurcation diagrams 
as a function of At with numerical basins of attrac- 
tion superimposed on a constant v- or w-plane. The 
other is the numerical basins of attraction with the 
stable asymptotes superimposed on the phase plane 
(u, v) with selected values of At. 

To obtain a bifurcation diagram with numeri- 
cal basins of attraction superimposed on the CM-2 
or CM-5, the preselected domain of initial data on 
a constant v- or rt-plane and the preselected range 
of the At parameter are divided into 512 equal in- 
crements. For the bifurcation part of the compu- 
tations, with each initial datum and At, the dis- 
cretized equations are preiterated 5,000-9,000 steps 
before the next 6,000 iterations (more or less de- 
pending on the problem and scheme) are plotted. 
The preiterations are necessary in order for the so- 
lutions to settle to their asymptotic value. A high 
number of iterations are overlaid on the same plot 
in order to detect periodic orbits or invariant sets. 
The reader is reminded that with this method of 
computing the bifurcation diagrams, only the sta- 
ble branches are obtained. While computing the 
bifurcation diagrams it is possible to overlay basins 
of attraction for each value of At used. For the 
basins of attraction part of the computations with 
each value of At used, we keep track of where each 
initial datum asymptotically approaches and color 
code them (as a vertical strip) according to the 
individual asymptotes. While efforts were made 
to match color coding of adjacent strips on the 


bifurcation diagram, it was not always practical or 
possible. Care must, therefore, be taken when in- 
terpreting these overlays. See Secs. 6.1 and 6.2 for 
discussions. 

For the basins of attraction on the phase plane 
(w, v) with selected values of At and the stable 
asymptotes superimposed, the (u> v) domain is di- 
vided into 512 x 512 points of initial datum. With 
each initial datum and At, we preiterate the re- 
spective discretized equation 5,000-9,000 steps and 
plot the next 6,000 steps to produce the asymp- 
totes (fixed points of various order and limit cy- 
cles). Again, for the basins of attraction part of the 
computations, for each value of At used, we keep 
track of where each initial datum asymptotically 
approaches and color code them according to the 
individual asymptotes. Details of the techniques 
used for detection of asymptotes and basins of at- 
traction are given in the appendix of Sweby & Yee 
[1991]. Note that in all of the plots, if color printing 
is not available, the different shades of gray repre- 
sent different colors. 

Selected results for both representations of nu- 
merical basins of attraction are shown in Figs. 6.1- 
6.33. In the plots, r = At. The “r=aDt” label 
denotes a scaling parameter “a” (set to unity for 
calculations presented here) times the time step At. 
White dots and white curves on the basins of attrac- 
tion with bifurcation diagrams superimposed repre- 
sent the bifurcation curves. White dots and white 
closed curves on the basins of attraction with the 
numerical asymptotes superimposed represent the 
stable fixed points, stable periodic solutions, or sta- 
ble limit cycles. The black regions represent diver- 
gent solutions. 

Note that the preselected regions of At and the 
selected At for the phase diagrams for both rep- 
resentations of numerical basins of attraction were 
determined by examining a wide range of At. In 
most cases, we examined At from close to zero up 
to one million. What are shown in Figs. 6.1-6.33 
represent some of the At ranges that are most in- 
teresting. Note also that for all of these models, 
using At = 0.1 is approximately 1/10 (or less) the 
time step limit of standard explicit methods [Yee & 
Sweby, 1992]. 

The streaks on some of the plots are either 
due to the nonsettling of the solutions within the 
prescribed number of iterations or the existence of 
small isolated spurious asymptotes. Due to the high 
cost of computation, no further attempts were made 
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to refine their detailed behavior since our purpose 
was to show how, in general, the different numer- 
ical methods behave in the context of nonlinear 
dynamics. From our numerical studies, the mid- 
point implicit method (linearized or iterative meth- 
ods) behaves the same as, or very similarly to, the 
trapezoidal method for the three models. Thus, 
no figures will be shown of the midpoint implicit 
method. In general, the dynamics of implicit LMMs 
are very different from the dynamics of standard ex- 
plicit methods. See reference [Yee &: Sweby, 1992] 
for the dynamics of nine explicit methods. 

6 . 1 . Numerical results for the 

dissipative complex equation 

Figures 6.1-6.10 show selected results for the two 
representations of numerical basins of attraction for 
model (4.1) for e = 1. The exact solution for (4.1) 
with e = 1 is a stable limit cycle with unit radius 
centered at (0,0). The “exact” basin of attraction 
for the limit cycle is the entire (u, v) plane except 
the unstable fixed point (0,0). 

Figures 6. 1-6.3 compare the noniterative with 
the iterative procedures for solving the nonlinear 
algebraic equations using implicit Euler, trape- 
zoidal, and 3-level BDF methods. Figures 6.4-6.10 
show selected results of the stable numerical asymp- 
totes with basins of attraction superimposed us- 
ing four different At by the three implicit LMMs. 
The red regions are the numerical basins of attrac- 
tion for the stable limit cycle except in Fig. 6.5 
{At = 2.65), in Fig. 6.6 {At = 1.5, 2, 4), in Fig. 6.8 
(At = 1.5), in Fig. 6.9 (At = 2, 2.5, 4), and in 
Fig. 6.10 (At = 1.515, 2.5). The green regions 
shown in Figs. 6.1 and 6.3-6. 5 are the numerical 
basins of attraction for the stabilized fixed point 
(0,0). Note how the implicit method turns the un- 
stable fixed point (0, 0) of the ODE system into a 
stable one for At > 1. 

To aid in the understanding of some of the re- 
sults shown in Figs. 6.4-6.10, the following gives 
some explanation on how to interpret the basins 
of attraction diagrams with the stable numerical 
asymptotes superimposed. All of the selected time 
steps At shown in Figs. 6.4-6.10 are based on 
Figs. 6. 1-6.3 where the bifurcation diagrams with 
the basins of attraction are superimposed. These 
time steps were chosen to illustrate selected features 
of the different bifurcation phenomena on the (u, v) 
plane. 


For example, it is easier to understand Fig. 6.8 
(At — 1.5) using “trapezoidal 4- modified Newton” 
if we look at the fourth plot in Fig. 6.2. Figure 6.2 
shows that the original limit cycle bifurcates into 
a “period 2 type” limit cycle near At = 1.25. Fig- 
ure 6.2 shows distinctively that both rings share the 
same basin (only one distinct solid red basin of at- 
traction). The lack of a solid basin of attraction 
in the third plot of Fig. 6.8 is due to the coloring 
algorithm, which requires a repetition of the limit 
cycle within the time of integration in order to dis- 
tinguish basins of attraction. When this repetition 
is not present the resulting coloring gives a crude 
indication of trajectories. Additional preiteration 
steps (many more than 9,000) would likely alleviate 
the problem. Note that we need all of the trajecto- 
ries corresponding to the 512 x 512 initial data to 
settle to within the prescribed preiterations before 
a solid basin results. The fourth plot of Fig. 6.8 to- 
gether with Fig. 6.2 hints at a rapid period doubling 
transition to instability. 

Figure 6.9 for At = 2 illustrates failure of the 
coloring algorithm to detect the basin of attraction 
due to an insufficient number of preiterations. The 
result is again a crude indication of the trajectories. 
Even though not correctly colored, the non-black 
region gives the size of the basin. Also, the third 
plot of Fig. 6.3 gives a clear indication of the size of 
the corresponding numerical basin of attraction for 
At = 2. 

Figure 6.10 for At = 1.515 illustrates the sur- 
prising presence of an embedded region of instabil- 
ity within the basin of attraction of the limit cycle 
(see also Fig. 6.3). The lack of a distinct (single red 
colored) basin is again an artifact of the coloring al- 
gorithm. Additional preiteration steps would likely 
alleviate this. 

Figures 6.1, 6. 3-6. 6, and 6.8-6.10 illustrate 
the situation where unconditionally stable LMM 
schemes, such as the implicit Euler and 3-level BDF 
methods, can converge to a wrong solution if one 
picks the initial data inside the green region, even 
though this region contains valid physical initial 
data for the ODE. Thus, even though LMMs pre- 
served the same number of fixed points as the un- 
derlying ODE, these fixed points can change type 
and stability. This phenomenon is related to the 
“nonrobustness” of implicit methods sometimes 
experienced in CFD computations. In these types 
of computations where the initial data are not 
known, the highest probability of avoiding spurious 
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Fig. 6.1 Fig. 6.2 
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Fig- 6.3 Fig. 6.4 







Basins of Attraction Basins of Attraction 

Dissipative Complex Eqn., c = 1 Dissipative Complex Eqn., e = 1 

Trapezoidal + Newton (Max. Iter. = 15, Tol. = .00001) Trapezoidal + Mod.Newton(Max. Iter. = 15, Tol. = .00001) 
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Fig. 6.7 Fig. 6.8 
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Fig. 6.9 Fig. 6.10 
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asymptotes is achieved when a fraction of the al- 
lowable linearized stability limit of At is employed. 

These diagrams also illustrate the unreliability 
of trying to compute a true limit cycle with any 
sizable At. This should not be surprising since 
the scheme only gives an 0(At p ) (p is the order 
of the scheme) approximation to the solution tra- 
jectories. In addition, since the limit cycle is not 
a fixed point, we would expect inaccuracies to be 
introduced. However, inaccuracies are not easy to 
detect in practice, especially when a numerical so- 
lution produces the qualitative features expected. 
Overall, the trapezoidal method and the midpoint 
implicit method (figures not shown) give more accu- 
rate solutions for the limit cycle. Note that all four 
LMMs except the implicit Euler are second-order 
accurate. 

All of the four LMMs and four solution proce- 
dure combinations exhibit spurious stable and un- 
stable asymptotes except spurious steady states. It 
is fascinating to see the dramatic difference in 
shapes and sizes of numerical basins of attraction 
for the different methods and solution procedure 
combinations compared with the exact basin of at- 
traction. The evolution of the numerical basin of 
attraction as At changes is very traumatic for all 
four LMMs and for the same LMM with different 
solution procedures. For larger At, the noniterative 
(linearized) implicit Euler, “implicit Euler + New- 
ton”, and “straight Newton” give the same numer- 
ical basins of attraction. Simple iteration behaves 
similarly to typical explicit methods (in terms of 
stability and the size of numerical basin of attrac- 
tion) for all of the four studied implicit methods 
and for models (4.2) and (4.3) as well. The main 
advantage of the simple iteration procedure over 
nonLMM explicit methods is that spurious steady 
states cannot occur. 

Numerical experiments were also performed on 
the dissipative complex model (4.1) using the adap- 
tive time stepping strategies (5.20) and (5.21) based 
on local truncation error control. Our studies indi- 
cate that these two variable time step control im- 
plicit and explicit methods (5.20) and (5.21) can 
alleviate the spurious dynamics for most calcula- 
tions. However, the allowable time step, determined 
by (5.20) or (5.21), is too small for practical us- 
age, especially for the explicit method (5.21). For 
the implicit method (5.20), the allowable time step, 
determined by (5.20), is slightly larger than the 
explicit method (5.21), but it is still impractical 
to use. Numerical experiments also indicated that 


with the wrong combination of starting time step, 
initial data and tolerance value, spurious dynamics 
could occasionally be produced. 

6 . 2 . Numerical results for the 
predator — prey equation 

Selected results for the two representations of 
numerical basins of attraction for the predator- 
prey model (4.2) are shown in Figs. 6.11-6.22. 
Figures 6.11-6.13 compare the noniterative with the 
iterative procedures in solving the nonlinear alge- 
braic equations using implicit Euler, trapezoidal, 
and 3-level BDF methods. Figures 6.14-6.22 show 
selected results of the numerical asymptotes with 
basins of attraction superimposed using four differ- 
ent Ats by the three implicit LMMs. Here, except 
for Figs. 6.14, At = 0.4 and 0.415, the green re- 
gions represent the numerical basins of attraction 
for the stable spiral (2.1, 1.98) and red regions rep- 
resent the numerical basins of attraction for the sta- 
ble node (0, 0). The numerical basins of attraction 
in Figs. 6.17 and 6.18 with At = 0.1 appear to be 
the same as the exact basins of attraction of the DE 
(4.2). The numerical basin of attraction by some of 
the LMMs for the fixed point (0, 0) is larger than 
the corresponding exact basin of attraction for the 
DE (4.2). (See Figs. 6.15, and 6.16 when At = 0.1.) 
In this case, the numerical basins of attraction for 
the divergent solution (black region) is smaller than 
the true one. The dramatic difference in shapes 
and sizes of numerical basins of attraction for the 
different methods and solution procedure combina- 
tions compared with the exact basin of attraction is 
even more fascinating than the dissipative complex 
model. 

All four LMMs and solution procedure combi- 
nations, other than simple iteration, change the two 
saddle points into stable or unstable fixed points of 
other types as illustrated in Figs. 6.14-6.22. For the 
implicit Euler, the two fixed points (2.1, 1.98) and 
(0, 0) are unconditionally stable and the stabilized 
fixed points (1, 0) and (3, 0) (saddles for the origi- 
nal ODE) are almost unconditionally stable except 
for small At. This is most interesting in the sense 
that the numerical basins of attraction for the sta- 
ble exact fixed points Ue of the model (4.2) by the 
implicit Euler method were permanently altered for 
At near or larger than 3, as illustrated in Figs. 6.15- 
6.17. It would be easier to interpret the results in 
Fig. 6.11 if one interchanged the yellow and green 
colors for At > 1. Our studies also indicated that 
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Fig. 6.17 Fig. 6.18 
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for larger At, the linearized implicit Euler, “implicit 
Euler + Newton”, and “straight Newton” give the 
identical numerical basins of attraction. See also 
Yee & Sweby [1992] for some discussion and illus- 
tration. Although the trapezoidal method, 3-level 
BDF, and the midpoint implicit method did not 
turn the two saddle points (1, 0) and (3, 0) into 
stable fixed points of different type, they did turn 
the two saddle points into unstable fixed points of 
different types. 

The evolution of the numerical basin of attrac- 
tion (as At changes) is very dramatic by each of 
the four LMMs. Take for example the trapezoidal 
method (Figs. 6.18-6.20), where the scheme be- 
comes effectively unstable for large At. The size 
of the numerical basins of attraction for the sta- 
ble fixed points Ue shrink to almost nonexistence. 
T his phenomenon might be one of the contribut- 
ing factors to the unpopularity of the trapezoidal 
method in CFD. The basins are so fragmented and 
small for large At that they are beyond the accu- 
racy of the CM-2 to resolve and no further attempt 
was made. A better approach in computing these 
types of basins is to use interval arithmetic or the 
enclosure type method [Adams, 1990]. 

Our results for the predator-prey model indi- 
cate that linearized implicit methods are more ro- 
bust and are more efficient than the other three it- 
eration procedures and have a higher success rate of 
leading to physically correct steady states. In gen- 
eral, it seems that if one uses a At far below what 
the linearized stability limit predicts, one has a bet- 
ter chance of avoiding spurious dynamics. Other- 
wise, the knowledge of numerical basins of attrac- 
tion is vital in avoiding spurious dynamics when 
using a fixed time step that is larger than the stabil- 
ity limit of the standard explicit methods studied in 
Yee &: Sweby [1992]. Comparing the current results 
with those in Yee & Sweby [1992], the implication 
is that unconditionally stable implicit methods are, 
in general, safer to use and have larger numerical 
basins of attraction than explicit methods. How- 
ever, one cannot use too large a time step since the 
numerical basins of attraction can be so small that 
the initial data has to be very close to the exact 
solution for convergence. 

Numerical experiments performed on the vari- 
able time step control (5.20) and (5.21) for model 
(4.2) also indicate that, although variable time step 
controls are not foolproof, they might alleviate the 
spurious dynamics most of the time. One shortcom- 
ing is that in order to avoid spurious dynamics, the 


required size of At is impractical to use, especially 
for the explicit method (5.21). 

6.3. Numerical results for the 

perturbed Hamiltonian equation 

Selected results for the two representations of 
numerical basins of attraction for the perturbed 
Hamiltonian model (4.3) (viscous Burgers’ equa- 
tion with 3-point central in space) are shown in 
Figs. 6.23-6.33 for e = 0.1. Figures 6.23-6.25 com- 
pare the noniterative and iterative solution proce- 
dures for solving the nonlinear algebraic equations 
using implicit Euler, trapezoidal, and 3-level BDF 
methods. Figures 6.26-6.33 show selected results 
of the numerical asymptotes with basins of attrac- 
tion superimposed for four different Af’s. In all of 
Figs. 6.23-6.33, red regions represent the numerical 
basins of attraction for the stable spiral (1/3, 1/3). 
The numerical basins of attraction in Figs. 6.28, 
6.30, and 6.33 with At = 0.1 appear to be the same 
as the exact basins of attraction. The numerical 
basin of attraction for (1/3, 1/3) is larger than the 
corresponding exact basin of attraction for At = 1 
by the implicit Euler and the 3-level BDF meth- 
ods. See Figs. 6.26, 6.27, 6.31, and 6.32. Note also 
that the possibility of the numerical basin of at- 
traction being larger than the exact one does not 
always occur when the time step is the smallest. 
For larger At the linearized implicit Euler, “implicit 
Euler -I- Newton”, and “straight Newton” give the 
same numerical basins of attraction. A conclusion 
similar to that arrived at in Sec. 6.2 for nonitera- 
tive and iterative procedures applies here. Since the 
midpoint implicit method behaves similarly to the 
trapezoidal method and the 3-level BDF exhibits 
behavior similar to that discussed in the previous 
section, our discussion is restricted to the implicit 
Euler and trapezoidal methods. 

Implicit Euler Method: 

This is yet another interesting illustration of the use 
of an unconditionally stable implicit method where 
in practical computations, when the initial data are 
not known, the scheme has a higher chance of ob- 
taining a physically correct solution if one uses a 
At restriction slightly larger than that for the sta- 
bility limit of an explicit method. Figures 6.23 and 
6.26-6.28 show the two representations of numer- 
ical basins of attraction using the implicit Euler 
method. These figures show the generation of stable 
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Basins of Attraction 
Viscous Burgers Eqn., G = 0.1 
Linearized Implicit Euler, Central Difference in Space 
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Fig. 6 . 26 . ( Continued ) 


spurious asymptotes for At > 1. As A t increases 
further, the size of the same numerical basin de- 
creases and becomes fractal-like, and new nu- 
merical basins are generated. The numerical 
basin of attraction for (1/3, 1/3) was permanently 
altered for At near or larger than 10. In general, 
the linearized implicit Euler is more efficient and 
less likely to converge to a spurious asymptote than 
the other combinations of LMMs and solution 
procedures. 

Trapezoidal Method: 

Figures 6.24, 6.29, and 6.30 show the two repre- 
sentations of numerical basins of attraction using 
the trapezoidal method. In a manner similar to 
the linearized implicit Euler method, this scheme 
has a higher probability of obtaining a physically 
correct solution if one uses a At similar to that of 
an explicit method. The numerical basins of at- 
traction for (1/3, 1/3) computed by the linearized 
trapezoidal method are much larger than the cor- 
responding exact basin of attraction for At < 2. 


Their sizes are bigger than the ones generated by 
the implicit Euler method with the same At 
values. The same behavior exists for the 3-level 
BDF method. As At increases, the shrinkage of 
the numerical basins of attraction is more dramatic 
than for the other two LMMs. In most cases with 
larger At, the allowable initial data required to 
avoid spurious dynamics is impractical to use since 
the “safe” initial data has to be very close to the ex- 
act steady-state solution. Again, this phenomenon 
might be one of the contributing factor to the un- 
popularity of the trapezoidal method in CFD ap- 
plications. For large At, the linearized trapezoidal 
scheme becomes effectively unstable due to the frag- 
mentation of the numerical basins of attraction (see 
Yee & Sweby [1992] for the plots). Again, due to 
the high cost of double precision computations, no 
further attempts were made for large At. The com- 
putation of these basins requires an interval arith- 
metic or the enclosure [Adams, 1990] type of math- 
ematical operation before a more precise behavior 
can be revealed. 
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Basins of Attraction 

Viscous Burgers Eqn., E = 0.1 (Cent. Diff. in Space) 
3-Level BDF + Mod.Newton(Max. Iter. = 15, Tol. = .00001) 


At =0.1 


At = 1 


At =1.5 


At = 3 


Fig. 6.33 


7. General Discussion of Numerical 
Results 

Studies showed that all of the four implicit LMMs 
exhibit a drastic distortion but less shrinkage of 
the basin of attraction of the true solution than 
standard explicit methods studied in Yee &; Sweby 
[1992]. In some cases with smaller At (near the 
linearized stability limit of standard explicit meth- 
ods), the implicit LMMs exhibit enlargement of the 
basins of attraction of the true solution. Overall, 
the numerical basins of attraction of a noniterative 
implicit procedure mimics more closely the basins of 
attraction of the continuum than the studied itera- 
tive implicit procedures for the four implicit LMMs. 
In general, the numerical basins of attraction bear 
no resemblance to the exact basins of attraction. 


The size can increase or decrease depending on the 
time step. Also, the possible existence of the largest 
numerical basin of attraction that is larger than the 
exact one does not occur when the time step is the 
smallest. The dynamics of numerics of the implicit 
methods differ significantly from each other, and the 
different methods of solving the resulting nonlin- 
ear algebraic equations are very different from each 
other since different numerical methods and solu- 
tion procedures result in entirely different nonlinear 
discrete maps. Although unconditionally stable im- 
plicit methods allow theoretically large At, the nu- 
merical basins of attraction (allowable initial data) 
for large At sometimes are so fragmented and/or 
so small that the safe (or practical) choice of At is 
slightly larger or comparable to the stability limit 
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of standard explicit methods. In general, if one 
uses a At that is a fraction of the stability limit, one 
has a higher chance of convergence to the correct 
asymptote. Studies in Yee Sz Sweby [1992] for stan- 
dard explicit methods confirmed this phenomenon. 

One of the causes of the above behavior of im- 
plicit LMMs is the existence of stable and unstable 
spurious asymptotes other than steady states which 
have a similar detrimental (in terms of robustness) 
effect as explicit methods. Another cause of the 
observed behavior is due to the fact that an unsta- 
ble fixed point can become a stable fixed point and 
can change type, e.g., from a saddle to a stable or 
unstable node. One consequence of the stabiliza- 
tion of unstable fixed points is a distortion, shrink- 
age and/or segmentation of the resulting numerical 
basin of attraction of the true steady states. An- 
other consequence of this behavior is that the flow 
pattern can change topology as the discretized pa- 
rameter is varied. Thus, even though LMMs pre- 
serve the same number (but not the same types) of 
fixed points as the underlying DEs, the numerical 
basins of attraction of LMMs do not coincide with 
the exact basins of attraction of the DEs even for 
small At, Some of the dynamics of the LMMs ob- 
served in our study can be used to explain the root 
of why one cannot achieve the theoretical linearized 
stability limit of the typical implicit LMMs in prac- 
tice when solving strongly nonlinear DEs (e.g., in 
CFD). 

Additional general remarks can be made for the 
following comparisons, independent of the four im- 
plicit LMMs: 

Simple Iteration versus Other Studied Methods : 

The stability and numerical basins of attraction by 
simple iteration are similar to those of standard 
Runge-Kutta and other commonly used explicit 
methods. The only advantage in using the “im- 
plicit method + simple iteration” over nonLMM 
explicit methods is that spurious steady states 
cannot occur. 

Noniterative versus Iterative : 

If one uses an implicit LMM for the time-dependent 
approach to obtaining steady-state numerical solu- 
tions, the linearized (or noniterative) version of the 
implicit methods are more efficient and less likely 
to converge to a spurious asymptote or get trapped 
in a spurious limit cycle than the other three stud- 
ied iterative procedures. Overall, the noniterative 


implicit Euler scheme is more stable and less likely 
to converge to a spurious asymptote than the 
other combinations of LMMs and solution proce- 
dures. The phenomenon can explain more precisely 
a contributing factor to the popularity of the non- 
iterative implicit Euler in CFD applications. 

Straight Newton versus Other Studied Methods: 

Studies indicated that contrary to popular belief, 
the initial data using the straight Newton method 
may not have to be close to the exact solution for 
convergence. Straight Newton also exhibits stable 
and unstable spurious asymptotes. Initial data can 
be reasonably removed from the asymptotic values 
and still be in the basin of attraction. However, the 
basins can be fragmented even though the corre- 
sponding exact basins of attraction are single closed 
domains. The cause of nonconvergence may just as 
readily be due to the fact that the numerical basins 
of attraction are fragmented. In many cases, the 
results obtained are better than those obtained by 
the trapezoidal and 3-level BDF methods (regard- 
less of the three iterative procedures). If one uses 
a time step slightly bigger than the stability limit 
of standard explicit methods for the four LMMs, 
straight Newton can have similar or better perfor- 
mance. In fact, using a large At by the linearized 
implicit Euler method or the implicit Euler + New- 
ton procedure has the same chance of obtaining the 
correct steady state as the straight Newton method 
if the initial data are not known or arbitrary initial 
data is taken. 

Variable Time Step Control: 

Our studies showed that the variable time step con- 
trol method (5.20) can occasionally stabilize un- 
stable fixed points, depending on the initial data, 
starting time step, and the value of “toll”. One 
shortcoming is that the size of At needed to avoid 
spurious dynamics is impractical to use, especially 
for the explicit method (5.21). 
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